On the Zero-Bias Conductance Peak in the Tunneling Spectroscopy 
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A generalized method of image, incorporated with the non-equilibrium Keldysh-Green's function 
formalism, is employed to investigate the tunneling spectroscopy of hybrid systems in the config- 
uration of planar junction. In particular, tunneling spectroscopies of several hybrid systems that 
exhibit zero-bias conductance peaks (ZBCP) are examined. The well-known metal-d-wave super- 
conductor (ND) junction is first examined in detail. Both the evolution of the ZBCP versus doping 
and the splitting of the ZBCP in magnetic fields are computed in the framework of the slave-boson 
mean field theory. Further extension of our method to analyze other states shows that states with 
particle-hole pairing, such as d-density wave and graphene sheet, are all equivalent to a simple ID 
model, which at the same time also describes the polyacetylene. We provide the criteria for the 
emergence of ZBCP. In particular, broken reflection symmetry at the microscopic level is shown to 
be a necessary condition for ZBCP to occur. 

PACS numbers: 74.20.-z, 74.50.+r, 74.80.FP, 74.20.Mn 



I. INTRODUCTION 

Since the pioneering work of Giaever,0 the tunneling 
measurement has been a major experimental method for 
investigations into the electronic states of condensed mat- 
ter systems.! 2 ! In the simplest setup, a metal with known 
spectral property is made in contact with a material X, 
forming an NX junction so that the electronic states 
of X can be probed. For many years, despite the fact 
that many insights into the spectral properties of many 
states have been gained from the differential conductance 
(dl/dV) curves obtained from tunneling measurements, 
nonetheless, unlike many other experiments, it is fair to 
say that there is no clear and solid statement as to ex- 
actly what bulk properties are being probed in tunneling 
measurements. For example, it is known that in neutron 
scattering experiments, the neutron intensity is a mea- 
sure of the imaginary part of the bulk spin susceptibility, 
Imx(k,w); no similar statement has ever been firmly es- 
tablished for tunneling measurements. 

The difficulty for establishing the relation between the 
tunneling conductance and the bulk quantities can be 
traced back to the very existence of the junction inter- 
face. It has been realized that the presence of the in- 
terface can change the conductance curve dramatically. 
A well-known example is the zero-bias conductance peak 
(ZBCP) observed in the tunneling spectra when X is a 
wave superconductor (ND junction) in (110) direction 
The appearance of the ZBCP is entirely tied up with the 
presence of the interface and its orientations, and there- 
fore can not be obtained by simple calculations based on 
bulk density of states. 

Recent theoretical analyses of the ZBCP have been 
mostly concentrated on the ND junctions. Furthermore, 
they are based largely on the standard BTK theoryp In 
the continuum limit, analytic expressions of the differen- 
tial conductance for general orientations of the interface 
were obtained. Numerical calculations were later car- 



ried out for the BTK theory in the lattice version! 
While these works have supplied insights into the ZBCP, 
they are, however, specifically designed for studying the 
ND junction. Moreover, because the relation of the con- 
ductance curve to the bulk quantities was not clearly 
manifested, essentially the numerical computation had 
to be done individually for each interface orientation. 
Another technical inconvenience is that the BTK theory 
is a mean-field theory based on solving the mean-field 
quasi-particle wavefunctions, it is thus difficult in this 
formulation to take into account the effects of interac- 
tion systematically. To extend into the study of other 
systems, especially those with strong correlations where 
almost all relevant models are on discrete lattices, it is 
therefore an urgent need to have a formulation which can 
go beyond the mean-field BTK formulations. As an il- 
lustration of our approach, in this paper we will focus on 
mean-field analysis of several tunneling problems. The 
effects of fluctuations and interactions will be discussed 
elsewhere. 

In this paper, we shall adopt an approach that is based 
on the non-equilibrium Keldysh-Green's function formal- 
ism. In the lowest order approximation, we will be able 
to express the differential conductance entirely in terms 
of bulk Green's functions and include the interface ef- 
fects. Thus, the relation of the conductance curve to 
the bulk quantities is clearly manifested. The tunneling 
between N and X will be treated as a perturbation, so 
that in the zeroth order the Green's function is the mean- 
field half-space Green's function that resides only on the 
semi-infinite plane and satisfies the boundary conditions 
to be specified later. Based on the half-space Green's 
function, <?, higher^xder corrections can be systemat- 
ically constructedEoB In particular, a class of infinite 
series in g, which consists of all elastic tunneling pro- 
cesses in the perturbation theory, will be considered and 
summed to alj-opiers for calculating the current across 
the junctiomQliS'til To fully take into account the tight- 
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binding nature of the problem, we shall employ discrete 
models for both the materials N and X and the tunnel 
junctions. Thus the essential quantity to be calculated 
is the half-space lattice Green's function for the X state. 
In resemblance to the conventional method of image, we 
express the half-space Green's function in terms of the 
bulk Green's functions propagating from the real source 
and a fictitious image source 
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(1) 



with the factor a accounting for the boundary conditions. 
In this picture the half-space Green's function is decom- 
posed into two parts: the real-source part comes solely 
from the bulk and hence reveals purely the bulk proper- 
ties, the image part contains all interface effects which 
are encoded in the factor a. In this way, the interface 
effects are clearly identified in the course of the analy- 
sis and one can pinpoint any departure from the bulk 
property. 

The factor a can be expressed in terms of bulk Green's 
function. Right on the interface, it is found 



a = G" 1 (d)G(-d). 



(2) 



Here d is an effective lattice constant whose precise mean- 
ing will be explained in below. Clearly, the tunneling 
spectrum can be classified according to whether the re- 
flection symmetry is broken or not. In the case when re- 
flection symmetry is broken with respect to the interface, 
one has G(—d) ^ G(d), hence c*o is not unity, possible 
zero modes may arise due to the presence of zeros in the 
denominator of the left hand side. The number of local- 
ized zero mode is thus determined by the order of zeros 
in the bulk Green's function G{d). In the lowest order 
approximation, the differential conductance is given by 
the local density of state at the interface 



dl/dV cx -J2lm{g (k,eV)}, 



(3) 



where go is g of Eq. (|l]) evaluated at the interface and e is 
the charge of an electron. Since cvo can be expressed en- 
tirely in terms of bulk Green's functions, this is then the 
relation between the bulk quantities and the differential 
conductance alluded to earlier. 

This paper is organized as follows. In Sec. |J, we out- 
line the theoretical formulation and derive the general- 
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ized method of image for discrete lattices. In Sec 
this method is applied to the study of tunneling spec- 
troscopies for various systems. We first study the ND 
junctions at various surface orientations and examine the 
doping dependence of the ZBCP using mean-field slave 
boson theory. We then study the effects of applied mag- 
netic fields perpendicular to the ab plane. A one dimen- 
sion model, based on the structure of polyacetylene, is 
then studied in Sec. 1IIC. On the basis of this model, we 



further apply this method to investigate tunneling into 
d-density-wave states and graphite sheets. We conclude 
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FIG. 1: (a) A typical configuration for the tunnel junctions 
studied in this paper: a (100) lattice on the left side connected 
to a (110) lattice on the right side. The dashed lines between 
the two lattices indicate hopping due to the tunneling Hamil- 
tonian. The effective one dimensional lattices obtained from 
Fourier transformation along the direction parallel to the in- 
terface are shown in below. In (b) and (c) we show explicitly 
the hard walls of semi-infinite lattices at (110) and (210) ori- 
entations, a and b indicate the crystalline axes. The boundary 
at x = cuts off the lattice and hence there can be no hopping 
across the boundary without the tunneling Hamiltonian. We 
implement this boundary condition by setting up hard walls 
at all lattice planes reachable by the boundary sites (filled cir- 
cles in the 2D lattices) and requiring the Green's function to 
vanish over these planes. In the (110) case, if only the near- 
est neighbor (n.n) hopping is considered then only the first 
hard wall is needed; the second hard wall imposes additional 
boundary conditions when there are next nearest neighbor 
(n.n.n) hopping. For (210) orientation, however, there are 
two hard walls even with only n.n. hopping. A third hard 
wall would be needed if one considers n.n.n. hopping in the 
(210) orientation. 



in Sec. [w] with some comments on the significance and 
further applications of our formulation. The Appendix 
describes techniques for deriving the current expressions 
for the tunnel junctions studied in the text. 



II. THEORETICAL FORMULATION AND 
GENERALIZED METHOD OF IMAGE 

A. Theoretical model 

We start by modeling the planar junction. As illus- 
trated in Fig. 0, the tunnel junction consists of two trim- 
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cated two dimensional (2D) lattices connected through 
a tunnel barrier, with the left half the normal (N) elec- 
trode (-co < x < —a^, ah is the lattice constant) and 
the right the test (X) electrode (0 < x < oo).E3 We take 
the interface the y direction. The total Hamiltonian of 
the system thus comprises two parts: the Hamiltonian 
Hq = Hl + Hr for the left and right electrodes, and the 
tunneling Hamiltonian which connects the surface points 
at x = —ah and x = 

Ht= J2 t (\Vl-Vr\)c\ a c rcT + ^. (4) 

Here a are spin indices, and yi, y r are the y-coordinates of 
the surface sites on the left and right electrodes; c/ CT , c ra 
are the corresponding electron annihilation operators, t 
is the tunneling amplitude whose magnitude models the 
barrier height in the tunnel junction. Since all points 
over the interface layers contribute to the tunneling pro- 
cess, one has to sum over all interface sites. Suppose the 
chemical potentials on the left and the right electrodes 
are and /ir, respectively, the total grand Hamiltonian 
is then given by 

K = (H L - h l N l ) + {Hr- m N R )+H T 

= K + H T . (5) 

The difference (/xl — [ir) is fixed to be the voltage drop 
eV across the junction. 

To calculate the tunneling current, we choose the un- 
perturbed state to be the ground state of Kq and adia- 
batically turn on Ht- In the Heisenberg picture, the tun- 
neling current is obtained from the time rate pf change 
of the particle number Nl of the left electrodetS (we set 
H = 1 throughout) 

I(t) = +ie([N L ,H T ] ) 

= +ie ]T {t(clc ra )-t*(clc la )}. (6) 

The expectation values (• • • ) here represent the ensemble 
average Tr[Z _1 exp (—(3Ko) ■■■]. 

In actual experiments, the normal metal on the left 
electrode could be in any orientations, and the detail 
connection between the two lattices may also cause com- 
plications in the tunneling spectroscopy. To be definite, 
however, in our model we fix the lattice on the left side at 
(100) orientation and connect its boundary sites to those 
of the right at x — (Fig. [j](a)). As one can observe 
easily, the system is translational invariant along the in- 
terface direction with period a^, the lattice constant of 
the left side. We exploit this symmetry by making a par- 
tial Fourier transformation along the interface direction 
in Eq. (^|) and arriving at 

I(t) = +ieJ2t(k y )(cl(k y )c ra (k v ))+h.c. (7) 

k y ,a 

The range of k y is determined by the periodicity of the 
interface sites along y direction, hence 

- Tr/a L < k y < ir/a L . (8) 



We emphasize that the problem is now effectively one 
dimensional: in Eq. (R) different k y modes are decou- 
pled completely. Moreover, only surface quantities are 
involved. These are very appealing features especially 
for the feasibility of our method of image, as we will dis- 
cuss in the following section. 

In the Keldysh Green's function formulation the time 
evolution of the density matrix can be. formally solved 
as a closed time-ordered path integraljlil the expectation 
value {c\ a {k y )c ra {ky)) in Eq. ([?]) is then related to the 
components of Keldysh's Green's functions over the close 
time-path. One can then calculate perturbatively the av- 
erage current I in terms of the zeroth order Green's func- 
tion. Details of this calculation can be found in Ref. |?] 
and an outline is presented in the Appendix. Here an es- 
sential difference from earlier work is that previously the 
Green's functions were obtained through directly solv- 
ing the equation of motion, while here we shall make use 
of the method of image elucidated in the following sec- 
tion. In this way the current approach is more general 
and versatile, and can be easily applied to various hybrid 
systems. 



B. Generalized method of image 

In our scheme for the calculation of the tunneling cur- 
rent, the building blocks are the zeroth order half-space 
Green's functions (see the Appendix). Because in the 
zeroth order, lattices on the left and right sides are dis- 
connected, the Green's functions are defined only for each 
semi-infinite plane. Therefore, lattice points on the in- 
terface will have "dangling bonds" . Effectively, as shown 
m Figs. §6) and (c), we are imposing hard-wall bound- 
ary conditions at the end points of these dangling bonds. 
One thus envisages a method of image similar to that in 
electrostatics. 

In the usual practice, the method of image is done 
for the continuum differential equations. It is based on 
the principle of superposition and the uniqueness of the 
solutions.li3 When applying it to the discrete lattice, one 
encounters the difficulty that the image point to any 
source point r, may not locate right at the allowed lat- 
tice points. To overcome this difficulty, we note that for 
each semi-infinite lattice there is discrete translational 
invariance along the surface direction, which we choose 
as the y direction. Furthermore, since in the analysis 
of tunneling problems each electrode is considered to 
be in steady states, time-translational invariance is pre- 
served in the individual half-space. In the subsequent 
sections we shall fully exploit these symmetries and thus 
will be concerned with the half-space Green's function 
in its Fourier space representation g(x, ), which 

effectively propagates from x' to x. For each k y and to 
one is therefore dealing with an effective one dimensional 
(ID) system (Fig. 0). 

As a demonstration of the method, let us consider a 
2D semi-infinite square lattice with lattice constant a ex- 
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tending over the region x > at orientation (hkO). The 
hard-wall boundary condition prescribes the half-space 
Green's functions to vanish over the hard walls, which 
consist of all points where the boundary sites can reach 
away from the bulk lattice (Figs. |l|(fc), (c)). For gen- 
eral surface orientations (hkO) and with only nearest- 
neighbor (n.n.) hopping one can find that the number 
of hard walls is given by max{|/i|, \ k\}. Let us consider 
first the single hard- wall configurations, which includes 
the (100) and the (110) orientations (when there is no 
next n.n hopping in the latter). As we shall discuss later, 
the multi-hard-wall problem are simple generalizations to 
the single hard- wall cases. 

For single hard-wall case the only hard wall is located 
at x = —d, where d = af\Jh 2 + k 2 is the spacing between 
two consecutive (hkO) planes. Since the Green's function 
must vanish on the hard wall regardless the position of 
the source point x' , one imposes the boundary condition 



g(-d,x';k y ;u) = 0. 



(9) 



To implement the method of image, we construct the 
half-space Green's function g(x, x'; k y \ u>) from the full- 
space Green's function G(x, x 1 ; k y \uS) as 

g(x, x'; k y ; uS) = G(x,x ; k y ;u>) 

— G(x,x 1 ;k y ;u>)a(x ;k y ;u>), (10) 

where x[ = —2d— a;' is the image point of the point source 
x' with respect to the hard wall x = — d. The Green's 
function G(x,x') describes direct propagation from the 
point source to the point x, while G(x,Xi) propagates 
from the image point to x. The factor is determined by 
fitting the boundary condition (||), we find 



a(x') = G -1 (-d,-2d- x')G{-d,x') 
= G- 1 {d + x')G{~d-x'). 



(11) 



Here and in the following we suppress the k y and u de- 
pendence whenever no confusion would arise. In going 
from the first to the second expressions above, we have 
used G(x,x') = G(x — x'), namely that the full-space 
Green's functions are translational invariant along the x 
direction. However, this is not essential for establishing 
the method of image. It is used here only for brevity. 
For systems without translational symmetry along x di- 
rection (such as the d-density-wave state to be discussed 



in Sec. HID), the following discussion still proceeds with 
only minor modification. 

From the half-space Green's function g(x,x'), one ob- 
tains the surface Green's function go by setting x = x' = 
0. In the Fourier space, go can be expressed by 

3o = G(k x ) x [1 - exp(2ik x d)a ] . (12) 

—ir/d<k x <lt/d 

Here ao = ct(0) does not depends on k x and the sum over 
k x extends over the first Brillouin zone of the effective ID 
lattice. The advantage of this formulation is clearly seen 
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FIG. 2: The method of image applied to d-wave superconduc- 
tors: the propagation (a) from the source A to the point B 
through the reflected path AOB in the presence of a hard-wall 
boundary can be replaced by (6) a direct path A'B emanating 
from a fictitious source at A' where the boundary is absent. 



from (jlj): the surface Green's function is obtained from 
combinations of full-space Green's functions. For differ- 
ent surface orientations, one simply rotates the full-space 
Green's function to the appropriate angle. Furthermore, 
it is also clear that here we have a scheme for studying 
the effects of interactions and fluctuations in tunneling 
problems. Essentially one can take these effects into ac- 
count through the bulk Green's function. Here, however, 
we shall concentrate on mean-field treatments and defer 
correlation effects to a separate publication. 

It is when dealing with lattices with an anisotropic or- 
der parameter that one could most easily appreciate the 
power of the present formulation. For instance in deal- 
ing with d-wave superconductors, apart from fitting the 
boundary conditions ([})), a also takes care of the different 
gap structures for propagation along the reflected path 
and the fictitious path (such as AO and A'O depicted in 
Fig. |). In the presence of reflection symmetry (such as 
an s-wave superconductor, or a d-wave superconductor 
at (100) orientation), since the gap structure as seen by 
these two paths are identical, the full-space bulk Green's 
function possesses the symmetry G(d) = G{—d). There- 
fore a becomes (independent of k y and to) universally 
equal to the identity matrix and Eq. (jl^) reduces to the 
familiar formQO 



3o = ^2 G ( fca; ) x 2sin 2 (fc x d) . 

— 7r/d<fe x <7r/d 



(13) 



For general orientations or when taking into account 
next nearest neighbor hopping, as noted earlier, there 
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could be more than one hard walls. In these circum- 
stances the surface Green's function must satisfy the 
boundary condition that it vanishes on all these hard 
walls simultaneously. This is a simple generalization of 
the single hard-wall problem. For instance, let us con- 
sider the (210) case with n.n hopping: as depicted in 
Fig. |l](c) there are two hard walls located at x = —d and 
—2d, where d = a/yS. Analogous to the single hard- wall 
problem, we write the half-space Green's function 

g(x,x') = G{x,x')-G{x,x' 1 )ai{x')-G{x,x' 2 )a 2 {x'){lA) 

with x\ = —2d — x', x' 2 = —id — x' being the location of 
the image sources, and a%, a 2 determined by the bound- 
ary conditions 



g(-d,x')=Q = g(-2d,x'). 



(15) 



In other words, for the point source at x' , each hard 
wall "generates" an image source on the other side of the 
surface and introduces an a factor which accounts for the 
additional boundary conditions. The half-space Green's 
function is a superposition of contributions from the real 
and all image sources. To obtain the surface Green's 
function, one again substitutes x = x' — into Eq. (|T^). 

Before proceeding to the applications in the following 
sections, we comment that the pr esent method is not re- 
stricted to square lattices. In Sec. Ill E we will apply this 
method to systems involving honeycomb lattices. Indeed 
our generalized method of image relies only on the possi- 
bility of reducing 2D lattices into ID structures through 
a Fourier transformation in the transverse direction. 



III. TUNNELING SPECTROSCOPY IN 
HYBRID SYSTEMS 

A. Normal metal— d-wave Superconductors 

We first study the a6-plane tunneling between a normal 
metal and a d-wave superconductor. The superconductor 
occupies the half-space x > and is described by the 
mean-field Hamiltonian 



Hr, 



+ ^2 A ^( c ^ c 3i ~ c a c it) + h - c - > ( 16 ) 
(v) 

where (ij) denotes the nearest- neighbor (n.n.) bond, 
(ij)' the next nearest neighbor (n.n.n.) bond, tR and t' R 
are hopping amplitudes between n.n and n.n.n sites, re- 
spectively; Ay is the mean-field pairing amplitude which 
possesses the d-w&ve symmetry 



A, 



A for r, = n + a, 
— Ao for Yj = Yi + b. 



(17) 



The normal metal on the left is modeled by a Hamiltonian 
similar to Hr but with only n.n hopping terms. 



To obtain the corresponding ID structure, we Fourier 
transform the Hamiltonian along the y direction. For ex- 
ample, at (110) orientation if including only n.n. hopping 
Hr becomes 



H T - 



E 



- /,I "" S ' ^f=) 4a( k y) c i+i<?( k y) 



k,. 



^ 2iA sin[ —j= ) [c^(k y )c l+n (-k y ) 



+Ci[(-ky)c i+ i^(ky)} + h.c. (18) 

Here Cj are the electron annihilation operators for the 
ID lattice at the z-th site (see Fig. [!](&)) and a is the 
lattice constant of the original 2D lattice. The lattice 
constant of the ID lattice is identical to the distance d 
between two consecutive (hkO) planes (cf. Fig. |l|). This 
ID structure of the problem is very helpful for us since 
each of the lattice planes (hkO) now becomes a point on 
the x axis. This enables us to define for each lattice site 
the corresponding image sites with respect to the hard 
walls£3 

It is convenient to use the Nambu notation which dis- 
tinguishes particle apd hole components. We define the 
spinor field-operatortH 



(19) 



The upper and the lower components of 'J',; correspond 
to the particle and hole components, respectively. For 
(110) orientation one can then write Hr of ( |l8| ) in the 
form 



with 



H, 



i,i±l 



-2i fl cos(^f) ±2zA sin(if) 
±2iA sin(^f) -2t fl cos(^f) 



(21) 



In order to apply the Keldysh formulation to calculat- 
ing the tunneling current, as detailed in the Appendix, 
the basic quantity one shall need is the (bare) half-space 
retarded Green's function. According to our method of 
image this can be obtained from superposition of the full- 
space Green's functions. Therefore our remaining task is 
to find the full-space Green's function. For general inter- 
face orientations (hkO) all that we need is to rotate the 
full-space Green's function to the appropriate angle and 
then build up the half-space Green's function based on 
the recipe outlined in Sec. [n| 

To find the full-space Green's function, we go over to 
the momentum space and express the Hamiltonian (|l6| ) 
in Nambu's representation 

Hr = Y £ k c L c k<7 + E( AkC kT c -ki + A ic c -kiC kT ) 



k,<r 
k 



-kl 



£k A k 



"kt 
-kl 



(22) 
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Here the quasiparticle dispersion £k and the gap function 
Ak are given by 

6k = — 2ifl[cos(k-a) + cos(k-b)] 

— 4j/ R cos(k-a) cos(k-b), 
A k = -2A [cos(k-a) - cos(k-b)] . (23) 

For (hkO) orientation the lattice vectors a = 
a(cos8, — sin0), b = a(sin 9, cos 9), where a is the lat- 
tice constant and 9 is the angle between the a-axis and 
the re-direction (thus tan# = k/h). The Hanhltonian H R 
in the form ( p2f ) is readily diagonalized, the quasiparticle 
excitation energy is found to be ±£k = ±\/ e k + A k . 

The full-space retarded Green's function can be ob- 
tained from 

G( Xi , Xj )= G(k„k y ,Lj)xe ik ^-^ , (24) 

-ir/d<k x <7r/d 

where G(k x , k y , u>) = \uj + ir\ — H R (k x ,k y ,uj with H R 
the matrix in the second line of ( p2] ) and r\ an infinitesimal 



h 
h 
h 
Ia 



Here the range of the k y is given by (^) , 

A a = i/(2nM, aa -g r ,J) (31) 

are the spectral weight matrices for the electrode a = 
{L,R}, and f(u) is the Fermi function (at zero temper- 
ature it is simply the step function 0(— u>)). The indices 
1, 2 in the Green's functions and the spectral weight ma- 
trices refer respectively to the particle and the hole com- 
ponents in the Nambu representation, t = t(u>,k y ) is 
the tunneling amplitude between the two electrodes. It 
is remarkable that the expression for I 2 here generalizes 
that found in Ref. and is applicable to any interface 
orientation. For the special cases considered in Ref. [t], 
where the surface Green's functions are symmetric (for 
(100) orientation) or antisymmetric (for (110) orienta- 
tion), Eq. fl2S| ) reproduces previous results. From these 
formulas one can clearly identify the contributions from 
each channel in the tunneling process. In particular, Ii 
is the contribution from single particle tunneling and Ia 



positive number. The half-space bare Green's function g$ 
is then obtained from the method of image. In the tunnel- 
ing problem, the tunneling Hamiltonian brings in tunnel- 
ing events between the two sides of the tunnel junctions 
which "renormalize" the half-space Green's functions (see 
the Appendix). In the Keldysh formulation this is ex- 
pressed as a perturbation series which can be re-summed 
to all orders in the tunneling amplitude yielding the 
renormalized half-space Green's functions S3 With the as- 
sumption that the renormalized advanced and retarded 
half-space Green's functions satisfy 

[9 r a ^ = 9% a (25) 

(a, (5 = {L, R} labels the electrodes), we can express the 
tunneling current as 

I = h+h + h + U, (26) 

where 



the Andreev reflection (thus Ia depends on the particle 
and hole components of the spectral weight matrix Al). 

We now present some of our results. Fig. || shows the 
tunneling spectra for (110) and (210) orientations at the 
doping levels 5 = 0.08, 0.14, and 0.20. Here we study the 
doping dependence by resorting to the mean-held slave 
boson theory for the t-t'-J model. The electron opera- 
tors c and are then essentially the spinon operators 
and the Green's function for spinons as well. The holons 
condense so that (b) = \/S. The mean-held parameters 
tR, t' R , Ao, and the chemical potential j.i R for each doping 
are calculated self-consistently.Q It is obvious from Fig. y 
that the ZBCP is significantly reduced in the (210) ori- 
entation. Interestingly, for (110) orientation the ZBCP 
decreases upon increasing doping while for (210) case it 
grows and then falls with doping. Another interesting 
feature in the tunneling spectra is the subgap structures 
near ±2Ao in the (210) case. These may have originated 
from resonances due to broken surface pairs, resulting 



/oo 
dLot 2 [f{u - eV) - f(uj)}A Lai {u> - eV)A R , n (u;)\l + tg r RLA1 (cj)\ 2 , (27) 
-OO 

/OO 
dci 2 [/( W - eV) - f(u))A LA1 (0J ~ eV) Ke{A R , 12 {uo)[tg a LRai (oj){l + tg RLtll (u))}} , (28) 
-oo 

/OO 
Aujt^fiio - eV) - - eV)A R , 22 (u)\9RL,i2(u)\ 2 , (29) 

-OO 

fry 

/OO 
d^i 4 [/(^ - eV) - f(co + eV)]A L , u (u> - eV)A L , 22 (u + eV)\g r RR ^ 2 {iu)\ 2 . (30) 
-oo 
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FIG. 3: The total differential conductance for (110) and (210) 
interfaces at dopings S = 0.08 (solid lines), 0.14 (dotted lines), 
and 0.20 (dash lines). The weak link is modeled by the in- 
terface hopping t(u) = exp(— y^wo — |w|)/r). Here we use 
u>o = 11 Ao and T = Aq. 



from the dangling bonds in (210) orientations.lla 

The ZBCP originates from zero-energy surface states 
(or the midgap states) due to Andreev reflections. In 
our formulation these states arise from singularities in 
the image contributions which manifest as poles in the a 
factors. In the presence of a single hard wall, the poles 
are determined by the zeros of the following factor when 
j] = 



(3{k y ) = det[G{d;k y ,uj = 0)]. 



(32) 



This produces singular behavior in the Green's functions 
and results in the ZBCP. In the (100) case, since ao is 
simply the identity matrix the surface Green's function 
( |l3| ) is regular at u> = thus there is no ZBCP. 



B. Tunneling into current-carrying 
superconductors 

We now consider the afc-plane tunneling from normal 
metals to current-carrying superconductors. In experi- 
ments one applies magnetic field along the c-axis of the 
superconductor, so that a screening current is generated 
over the ab plane. When a quasiparticle tunnels across 
the surface layer, it acquires additional energy from the 
supercurrent. Thus the zero-energy surface state evolves 
in this case into two surface states with non-zero energy. 
In the tunneling spectra this appears as "splitting" of the 
ZBCP (Fig. |j). Fogelstrom et al. have analyzed splittings 
of the zero-energy peak in the surface density of states 
under applied fields in the continuum limit. tj Here we ex- 
amine the tunneling spectra base on our discrete model. 
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FIG. 4: Splitting of the ZBCP for various values of q y (for 
5 = 0.16). 



To marry with formulations in the previous section, we 
note that in the.presence of supercurrent the gap functi 
is modified asE3 



ion 



A, 



Aij exp (iq • (r^ 



,0), 



(33) 



where q = (0, q y ) is the superfluid momentum and is 
proportional to the magnetic field. We shall assume that 
tunneling events take place only within a shallow layer 
of order about the penetration depth from the surface, 
so that q y is approximately uniform in the region of our 
concern. This additional phase can be absorbed into the 
electron operator by the transformation ci„ — > a a exp(iq- 
rj). In Fourier space the Hamiltonian becomes 



H, 



k . cr 

E(4 



£k+q Ak 

kT u -kl )\ A* -e k _ 



^(Ak4 T c -ki + A kC-kic kT ) 

(34) 



'-k-l 



Here the quasiparticle dispersion e k aud the gap function 
Ak are those of Eq. ( p3| ) . After diagonalizing Hr above 
one finds the quasiparticle excitation energy becomes 



(±) _ 




A k .(35) 



The momentum space Green's function that is fed into 
Eq. ( p"2| ) is obtained in the same way: G(k x ,k y ,Lu) — 
[u> + ir/ — H R,{kx,k y , , with Hr the matrix in the 
second line of (|34|). 

Fig. ^ shows typical tunneling spectra for the split- 
ting of the ZBCP when increasing q y . Note that the 
slightly asymmetric splitting originates from the particle- 
hole asymmetry in e k . Fig. ^| plots the magnitude of the 
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FIG. 5: The dependence of splitting on magnetic field for dop- 
ing S = 0.12. The empty and full symbols represent data cal- 
culated, respectively, with and without self-consistently tak- 
ing into account the magnetic fields in solving the t-t'-J slave 
boson mean-field equations. In the former case, the super- 
conducting gap is strongly suppressed when q y > 0.65, where 
difficulty in convergence of the mean-field solution arises. 



splitting versus the applied magnetic field for underdoped 
case. For small q , the expansion of Eq ( |35|) leads to linear 

as observed in small applied 



splitting in the lowest order terms: E^ 
where E^ 



= ^ 



fields. For higher fields, one has to retain the full q 
dependence, resulting in the bending of the splitting. 
This is purely due to the lattice effect. Also shown in 
Fig. U are the results taking into account suppression 
of the superconducting gap under magnetic fields self- 
consistently. The curve is seen to be "pushed" inwards 
while maintaining similar features. Note that quantita- 



tive agreement with experimental observationsEinij can 
be obtained by fitting scales of our results to the experi- 
mental data. Nevertheless, we did not observe any zero- 
field splitting at overdoping. This is in contrast with 
the experiment of Ref. |2^. The mechanism inducing this 
splitting might have eluded from our simple model. 

The doping dependence of splitting is also shown in 
Fig. ^ for q y = 0.15 and 0.80, which are respectively in 
the linear and the saturated regimes in Fig. |. Note that 
the splitting increases with doping, in agreement with 
Ref. || 

In passing we point out that the splitting depends sen- 
sitively on the Fermi surface topology. Indeed for /ir = 
we find no splitting of the ZBCP whatever the value 
of q y is. One can confirm this analytically by making 
an asymptotic expansion of the Green's function around 
uj = 0. At fin — one finds the conductance peak invari- 
antly stays at the zero bias. 
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FIG. 6: Splitting versus dopings for q y — 0.15 (open squares) 
and 0.8 (solid squares). 
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FIG. 7: The structure of polyacetylene and the ID ti-t2 
model. Filled and empty circles are lattice points over the 
A and B sublattices. 



C. Polyacetylene 

Up to this point, we have considered tunnel junctions 
with superconducting test electrodes, where particles of 
opposite spins form pairs. In this and the following sec- 
tions we will consider systems which exhibits particle- 
hole pairing over bipartite lattices. To start with we shall 
consider first a simple ID model based on the structure 
of polyacetylene .E3 This will turn out to be very helpful 
for understanding results in the following sections. Most 
importantly, it provides the criteria for the formation of 
midgap states in semi-infinite bipartite systems. 

The model we shall examine here is a ID chain with 
alternating hopping amplitudes ti, ti as shown in Fig. [?|. 
The separation between the lattice points is taken to be 



a constant ak-i It is convenient to categorize the lattice 
points into A and B sublattices and express the Hamil- 
tonian for this "ti-t2 model" as 



H, 



E 



+ At 



t 2 c 



-i+l 



h.C. 



(36) 
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Here cf annihilates electrons over site i on the a = 
{A, B} sublattice (spin indices a will be omitted through- 
out), and the sum run over sites i in the B sublattice only. 
Going over to momentum space, one finds 



Hr 



^A^cjf + h.c. 



( ^ cf ) 



A fc 




(37) 



where 



A, 



-t\ exp(ifca) — ti exp(— i/ca) 

-{ti + ti) cos(fca) — i(t\ — ti) sin(fca) . (38) 



The Hamiltonian ([37|) can be diagonalized easily and 
the quasiparticle excitation energy is found to be ±|Afc|. 
Note that the real part of A& is similar to the usual hop- 
ping energy in one dimensional chain. Therefore, when 
Im{Afc} oc (tx — tz) 7^ a single particle excitation gap 
opens at the chemical potential. 

For semi-infinite chain, there are two possible config- 
urations with the terminating site being an A or a B 
sublattice point. In either case we choose the boundary 
point the origin x = and construct the surface Green's 
function utilizing the method of image 



So 



G(0, 0) - G(0, -2a)G" 1 (-a, -2a)G(-a, 0) . (39) 



In this formula, the appropriate Green's functions should 
be used depending the the type of the end point; for 
instance, in the case of an A- type boundary even/odd 
sites are attributed to the A/B sublattices. The retarded 
Green's function for infinite system is 

G(x?,t;x?,0) = -iQ(t) ({c?(t),cf (0)}) (40) 

for xf over the a sublattices; the braces denote the anti- 
commutators. After Fourier transformation in time, the 
Green's function is obtained from 



G(x?,x$;w)= J2 G af} (k,Lu)xe^-^. (41) 

-7r/2a<fc x <7r/2a 



Utilizing Hr the matrix in the second line of (|37|), we 
write for brevity the momentum space Green's functions 
in matrix form 



G a p{k,u) = [oj + irj- H R {k,Lu)}~p 

1 / oj + ir) A fe 

(w + if]) 1 -E 2 k \ A* k UJ + i r] 



(42) 



a/3 



The matrix elements are assigned to the Green's func- 
tions according to the convention in (|37|), namely a — 
A, B corresponds to a = 1, 2 respectively. 

Since in Eq. ( (39| ) the surface Green's function g is ex- 
pressed as a combination of bulk Green's functions, the 



only possible source of singular behavior in go resides in 
the inverse part G _1 (— a,— 2a). In other words, the ex- 
istence of the zero-energy mode depends on the behavior 
of G(— a, —2a) at lo = 0. This is analogous to the ND 
junctions where the ZBCP results from the zeros of the 
determinant fl(k y ), Eq. (^). 

For example, in the case of A-tyj>e boundaries we find 
explicitly (setting rj = 0) 



Gba(—o>, -2a;cj = 0) 



a 

2^ 



7r/2o 



dfc- 



1 



7T/2Q h + h exp(-2ifca 

if ti < t 2 , 
l/2ti if *2 < h. 



(43) 



Here the index U BA" denotes the same meaning as in ( |42| ) 
and is used for emphasizing the correct Green's function 
to be used. From (43), when t\ < ti a sharp singularity 
in the surface Green's function go arises at uj — due to 
the divergent factor G^ A in (|3^) . On the other hand, for 
t\ > ti since Gg A is finite at u> — 0, no singular behavior 
in go could occur there. Thus one expects ZBCP in the 
former case while none in the latter. In the following 
sections we shall see that this provides for 2D bipartite 
systems a criterion for the range of transverse momenta 
where zero-energy states exist. For £>-type boundary the 
analysis is identical, except an exchange in the roles of t\ 
and ti. Therefore when the ZBCP shows up in an A- type 
chain, it must be absent in a B-type chain, and vice versa. 
This is shown in Fig. || for the case of polyacetylene. The 
current expression here is identical to Eq. (|l9|) given in 
the following section, except the extra sum over k y there. 



D. Normal metal— d-density wave states 

In underdoped cuprate superconductors, it is observed 
in experiments that there are signatures of a "par- 
tial" gap well above the superconducting temperature 
T c . This anomalous regime in the phase digram of the 
cuprate, superconductors is thus termed the pseudogap 
phase£3 Experiments also find that the pseudogap is con- 
sistent with a d-wave structure. Recently Chakravarty et 
al. proposed that the pseudogap phase of the underdosed 
cuprate is possibly the d-density-wave (DDW) stateB3 It 
is therefore of interest to examine the tunneling spectra 
of normal-metal-d-density-wave (N-DDW) junctions. 

The DDW state is characterized by the staggered flux 
in the elementary plaquettes of the lattice. The bond 
currents circulating the unit cell of the underlying square 
lattice break, among other symmetries, the invariance of 
translation by one lattice spacing and lead to a bipar- 
tite structure (Fig. ||). Obviously, if the interface cuts 
at (110) direction, the reflection symmetry is broken - 
in contrast to the (100) case. Therefore, we shall ex- 
amine the (110) direction with the following mean-field 
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direction in H Rl one finds 

H *= E { Hi-l4\(k y )cf(ky) 

0.10 - - xf,k y ,a 

+ A hl+1 cf f (k y )cf +1 (k y ) + h.c.} , (45) 
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FIG. 8: Typical tunneling conductance curves for polyacety- 
lene with A type (solid line) and B type (dashed line) end 
points. Here we take t\ = 2.25 eV and = 2.85 eV; 
thus the bandwidth is ti + %i = 5.1 eV and the gap width 
|tl — ta| = 0.6 eV. The linear chain on the left side has been 
taken a wideband material. In the tunneling Hamiltonian Ht 
we take t — 0.3. 
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(a) (b) 

FIG. 9: (a) The configuration of a 2D square lattice with d- 
density-wave order at (110) orientation and the corresponding 
ID model. Filled and empty circles label the A and B sub- 
lattices; the arrows indicate the directions of bond currents. 
Dashed lines extended from the boundary sites depict cou- 
pling to the left electrode through the tunneling Hamiltonian. 
(b) shows explicitly the bond variables in a doubled unit cell. 



HamiltonianEZI 

+X* (c£ b cf + cfj b cf) +h.c.} , (44) 

where cf annihilates an electron at site i over the a 
sublattice, and \ is the hopping amplitude on the bond 
(Fig. 0(6))- Making Fourier transformation along the y 



where A^±i = 2Re{xe ±tkya ^} with a the lattice con- 
stant of the square lattice; cf(k y ) is the electron annihi- 
lation operator for fixed k y at site xf over the ID lattice. 
Going over to momentum space, one finds, similar to ( |3^ ) 

with Ak = £k + «Ak- £k and Ak are given by Eq. (|23| ) 
with t R = -Re{x}, A = -Im{x}, and t' R = 0. The 
quasiparticle excitation energies is then obviously ±£k = 

±|A k | = ±y4TAj-. 

To find the tunneling current we apply again the 
Keldysh formulation outlined in the Appendix. For fixed 
k y the full-space retarded Green's function between the 
sites xf and Xj pertaining to the a and /3 sublattices is 
given by 

G(xf,t; x$, 0; k y ) = -ie(t)({cf(k y ,t), cf (k y , 0)}), (47) 
which in Fourier space for fixed frequency to becomes 
G(x?, X P)= J2 Gap(k x ,k y ,Lu) x e ikm( - x ?~ x iK (48) 

-7r/2d<fc a: <7r/2(i 

Here d = is the lattice spacing of the ID lattice and 

we have suppressed the u> and k y dependences on the left 
hand side. G a p has the same form given in Eq. ( [42] ) 
except that now a — A or B. 

From Ga/3, the half-space surface Green's function is 
obtained again using the method of image. The cur- 
rent expressions here, however, are distinct from those of 
(f27|)-(|30"|). Indeed since we are dealing with a single com- 
ponent Green's function the calculation is much simpler 
than previously. As shown in the Appendix, the current 
expression is here 

/oo 
&ut 2 [f{u-eV)- /(w)] 

xA L (u-eV)A R {u)\l + tg r RL {u)\ 2 . (49) 

This is exactly the single-particle current I\ of Eq. ( [27| ) 
for ND tunneling. There is no contribution from "An- 
dreev reflections" in N-DDW tunneling. This is due to 
the fact that in the DDW state the pairing takes place 
between particles and holes of momenta k and k + Q, 
with Q the nesting vector of 2D square lattices. Thus 
the Andreev reflected particles are still electrons whose 
response to the bias voltage are the same as the incident 
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FIG. 10: Typical conductance (dl/dV) curves for N-DDW 
junctions at (110) orientation in the (a) absence and (b) 
presence of in-plane magnetic field. Here the boundary sur- 
face consists of A sublattice sites and \ — (~tR ~ *Ao) = 
(-0.447 - O.li), rj = 0.01. Also shown in (6) are contribu- 
tions from the spin-up (dashed line) and spin-down (dotted 
line) components. The Zeeman splitting is here 0.24 Ao. The 
weak link is modeled by the same expression as in Fig. ^[ 

particles; as a result their contributions to the tunnel- 
ing current cancel exactly. In the ND junction, however, 
a particle is Andreev reflected as a hole, which behaves 
oppositely under applied bias. Fig. 0(a) shows a typ- 
ical plot for differential conductance versus voltage for 
N-DDW junctions. The conspicuous ZBCP agrees rwjth 
recent calculations done by Honerkamp and Sigrist.Ea 

The reason for the ZBCP here can be understood on 
the basis of the results in the previous section. Just like 
polyacetylene, the midgap states arises when go is singu- 
lar due to the zeros in the Green's function such as in 
Eq. (fi"3|). For each k y Eq. ( [i"5|) resembles the ti-t2 model 
with t\ = — Aj^-i and t% = — Aj.j+i. Therefore, for ex- 
ample, for A-type boundary one expects midgap states 
for the range of k y where 

A M _i > A i)i+1 or Im{x} sin (7^=) > ■ ( 50 ) 

Since here Im{x} = — Ao < 0, the above equation leads 
to -y/2n/a < k y < 0. 

We have so far considered only the case of vanishing 
chemical potential fiR in the DDW state. At finite chem- 
ical potential the grand Hamiltonian for the DDW is 
Kr = Hr — /irNr. Since the number operator 

Hence —/ArNr is diagonal and it simply shifts the ex- 
citation energy to E^ — (j,r. It is easy to check 



FIG. 11: Same as Fig. |lfj(a) but with next n.n hopping am- 
plitudes t' R = 0.0 (solid line), -0.03 (dotted line), and -0.06 
(dashed line). 



that this change in the excitation energy induces a shift 
u> — ► uj + fin in the Green's function. This is in sharp 
contrast with the ND case; there the chemical poten- 
tial shifts the quasiparticle energy ek — > ek — hr in the 
Green's function but not the frequency. This results in 
the distinct behavior of the ZBCP for ND and N-DDW 
junctions at finite [iR . 

For N-DDW junctions since u> — * cv+fj,R at finite chem- 
ical potential, the conductance peak is shifted from zero 
bias to the opposite value of the chemical potential —/jr. 
For ND junctions, however, the midgap state stays at 
lu = even at finite chemical potential, thus the con- 
ductance peak always position at zero bias (see Fig. ||). 
This shift has an obvious implication: the peak will split 
due to the Zeeman splitting (see Fig. |l^). The orbital 
effects of magnetic fields can be included by changing \ 
into ^e lq '( ri_r " for any nearest neighbor sites i, j. This 
takes into account the current induced near the interface. 
Since under this change both ek and Ak undergo shift- 
ing of k by q which can be absorbed into the summation 
of k, the peak does not split. Therefore, the splitting 
of ZBCP turn out the same for both in-plane and per- 
pendicular magnetic fields. This is in contrast to the ND 
junction where orbital effects dominate for perpendicular 
fields. 

In closing this section we note that since the next n.n 
term —At' R cos(k • a) cos(k • b) couples only lattice sites 
within each sublattice, it introduces only diagonal terms 
to (pff). Therefore, similar to the chemical potential, the 
next n.n terms cause the ZBCP and the spectrum to 
migrate when t' R 7^ 0. This is displayed in Fig. 11. 
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FIG. 12: Graphite sheets with (a) zigzag and (b) armchair 
boundaries and the corresponding ID models after k y trans- 
formation. Filled and empty circles represent respectively the 
A and B sublattices; a, b are lattice vectors. The dashed 
lines draw from the boundary sites indicate connections to 
the left electrode through the tunneling Hamiltonian similar 
to Fig. go). 



E. Graphite sheets 

So far we have considered systems involving only 
square lattices. As commented in the end of Sec. [nj our 
formulation is indeed quite general and can be applied to 
any systems which can be projected into ID structures. 
As an example, we study in this section the in-plane tun- 
neling from a normal metal into semi-infinite graphite 
sheets (NG junctions). 

In the 2D graphite sheets, each carbon atom is bonded 
to its three n.n via three electrons in the sp 2 orbital, 
forming a honeycomb lattice. The fourth electron (the 
7r electron), which resides in the p z orbital perpendicular 
to the 2D layer, is not bonded and is free to hop from 
site to site. In the tight-binding limit the Hamiltonian 
for the bulk graphite sheet is thus 

Hn =E -Toc^cf-Ttc^cf-^cf^.cf +h.c. (52) 

Here the lattice is divided into A and B sublattices, and 
a, b are the lattice vectors illustrated in Fig. [jj]. cf are 
electron annihilation operators at site i over the a sub- 
lattice, and 7^ are the hopping integrals. For simplicity 
we shall take 7o = 7i = 72 in the following. We will 
be interested in two orientations of the lattice: one with 
zigzag and the other with armchair boundaries. 

We first consider the zigzag case and choose the frame 
of coordinates as shown in Fig. |l2](a). Fourier transfor- 
mation in the transverse direction leads to ID Hamilto- 
nian which resembles ( |36| ) 

H R = 5>*ic£l (ky)cf (ky)-t 2 C^ {k y )cf +1 (ky) +h.C. 

(53) 



with here 

t\ = 270 cos(-^-fcj,et) and t 2 — 70 ■ (54) 

Further k x transformation brings Hr into the same form 
as (H|) with 

A k = -7o[e"^ a + 2cos(^fc y a)e^ a / 2 ] . (55) 

In applying the method of image, we note that the 
projected ID lattice for the zigzag case has alternating 
bond length, which breaks the reflection symmetry and 
hence implies the possible existence of the ZBCP. The 
alternating bond length, however, seems to cause diffi- 
culty in locating the image point of an arbitrary source 
site. For instance, the usual choice - the mirror image 
- does not always put the image point right on the lat- 
tice. Nevertheless, since in ID the hard wall becomes a 
point, as long as the Green's function propagating from 
the real source to the hard wall can be canceled by that 
from a fictitious source so that the boundary condition 
is satisfied, uniqueness of the half-space Green's function 
implies that the location of the fictitious source can be 
chosen at will. Indeed, this can be explicitly checked 
numerically. To be definite, we shall place the fictitious 
source at x — — (3/2)a and apply the method of image. 
The boundary condition <?(—&, x') — for all x' immedi- 
ately leads to 

90 = Gaa(0,0) 

- G A A(0,-3/2a)G^(-a,-3/2a)G^(-a,0). (56) 

Here we have labelled the attributes of the lattice points 
explicitly in the subscripts of the Green's functions. Just 
like polyacetylene, the midgap states arises when go is 
singular, namely at the zeros of Gba(~o-, — 3/2a) when 
■q — 0. From ( p4|) the correspondence to the t\-t2 model 
indicates that midgap states exist for k v which satisfy 
cos(fcj,\/3a/2) < 1/2, ie. when setting \/3a = 1 

- 7T < k y < -2vr/3 and 2tt/3 < k y < tt . (57) 

This is exactly what is found in band structure 
calculations EB 

For the zigzag orientation, apart from the zigzag 
boundary, there could also be the "bearded" boundary 
where the surface layer consists of B sites. This is remi- 
niscent of the case of -B-type end point of the t\-t2 model. 
Similar analysis as above can also carry over here with 
minor change. We find in this case the zero energy state 
arises from the range (v^a = 1) 

- 2tt/3 < k y < 2tt/3 . (58) 

The current expression are the same as Eq. (^9|) for N- 
DDW junctions. The corresponding tunneling spectra 
are solid and dashed lines shown in Fig. Ol 
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FIG. 13: Typical tunneling conductance curves for NG junc- 
tions with zigzag (solid line), bearded (dashed line), and arm- 
chair (dotted line) boundaries. Here 70 = 0.1 and the left 
electrode has been taken a wideband material. The weak 
link is modeled by the same expression as in Fig. ^ with here 
ujq = II70 and F = 70. 



Let us now consider the armchair case (Fig. 12(6)). 
After the Fourier transformation along the interface, one 
finds 



Hr = ^70 



c?\k y )cf(k y )e- ik * a 

+ cf\k y )ct 1 {k v y kya/2 +^-) 



(59) 

Note that here for each site Xi there are both A and 
B components as shown in Fig. |l^(6). The propagation 
between Xi and Xj thus compose of four components and 
the full-space Green's function form a 2 x 2 matrix 



(G(xf,xf) G{x?,xf) 



(60) 



Further k x transformation yields the Hamiltonian ([46] ) 
with here 



A k = -7o[e ifc » a +2cos(— k x a)e 



(61) 



Note that A(—k x ,k y ) — A(k x ,k y ) implies that the re- 
flection symmetry is not broken. Similar to the case of 
polyacetylene, in momentum space, G(k, lu) has exactly 
the same form as Eq. (fl2|). However, now the hard- wall 
boundary condition becomes a matrix equation 



— d 



for all a,0 = {A, B} , (62) 



where d — (\/3/2)a is the lattice constant of the projected 
ID lattice. The surface Green's function then takes the 
form 



.90 



G(0, 0) - G(0, -2d)G" 1 (-d, -2d)G(-d, 0) . (63) 



Since transnational symmetry is preserved in this ID lat- 
tice, we have G{-d, -2d) = G(d) and G(-d, 0) = G{-d). 
Reflection symmetry implies G(d) — G(—d) and conse- 
quently 



.go = G(0) - G(2d) 



(64) 



In this case, one thus expects no midgap states. 

Without loss of generality, we connect the A sublattice 
to the left side (Fig. 0(6)). The current expression is 
then the same as Eq. (|49[) , where go is replaced by the 
11 (or A A) component of the right hand side of Eq. (|64|). 
The dotted line of Fig. [l^ shows the conductance curve 
for the armchair case. 



IV. SUMMARY 

In summary, the generalized method of image that we 
developed has allowed us to deal with various tunneling 
problems in a unified manner, with full tight-binding na- 
ture being taken into account. In particular, we have ap- 
plied it in this paper to examine in-plane tunneling spec- 
tra of normal metal-G?-wave superconductor junctions, 
with and without external magnetic fields, at arbitrary 
crystalline orientations. The doping dependence of the 
ZBCPs is also studied within the mean-field slave bo- 
son approach. Our results for the splitting of ZBCP un- 
der applied magnetic field agrees well with recent experi- 
ments. We further showed that tunneling into G?-density- 
wave state at (110) orientation should display a sharp 
conductance peak at the chemical potential in the tunnel- 
ing spectra. This peak will shift away from the chemical 
potential if the next nearest neighbor hopping t' R exists, 
which also offers a way to measure t' R . Under in-plane 
magnetic fields, it also splits due to Zeeman splitting. 
These provide signatures to be looked for in experiments, 
especially in normal-metal-pseudogap-cuprate junctions 
for testing the proposal of Ref. |2(| As a demonstration 
of the general applicability of our formulation, we fur- 
ther consider tunneling into graphite sheets at the zigzag 
and armchair orientations. ZBCP is found in the zigzag 
case while no ZBCP should be displayed in the armchair 
case, coflsistent with findings in the study of graphite 
ribbons .123 We analyze these results on the basis of the 
ID t\-t2 model and obtain the criteria for the emergence 
of the ZBCP. 

The merit of our formulation lies in two aspects. 
Firstly, it offers a unified method for theoretical study 
of the tunneling spectroscopy. Recently, applying this 
method, we discover a remarkable even-odd effect in 
semi-infinite ID chains: for hopping amplitudes with 
even cycles there can be zero-energy localized edge-states, 
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while for those with odd cycles there is none@ Secondly, 
as already pointed out at the beginning of the paper, 
for the first time, our method allows us to express what 
is being measured in tunneling experiments in terms of 
bulk Green's function. For instance, in a single hard-wall 
configuration, suppose the normal metal on the left is a 
wideband material and the junction barrier is very high 
(t <C 1 in Eqs. (|7])-(|3(]) and (49)), tunneling experiments 
essentially measure the quantity 



dl/dV ex -^Im{<7o(k,eT/)} 

= - Im{G(k, eV) X [1 - exp(2i/c a; d)ao]} . 

(65) 



k CT 



Namely it probes the surface density of state, which is 
the bulk Green's function modulated by the factor in the 
square bracket. In our formulation the surface Green's 
function, and hence the surface density of states, decom- 
poses into two parts: one from the real source and the 
other from the image source. The real-source part re- 
sults purely from the bulk and hence reveals faithfully the 
bulk property; the image source part contains all surface 
effects and hence is responsible for any complications. 
In the presence of the reflection symmetry, we find that 
Q'o = 1 and the image part contribute another bulk term. 
Thus in this case the conductance curve exhibits only the 
bulk property (with, however, the van Hove singularity 
"rounded" by a sine factor as in Eq. (|l3|)). When the re- 
flection symmetry is broken, singular behavior may arise 
from the image part. This singularity is contained in 
ao and originates from the zeros of the Green's function 
(Eq. (]43|)) or its determinant when considering supercon- 
ductors (Eq. <M 
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APPENDIX A: CURRENT EXPRESSIONS 

In this appendix we outline techniques for calculating 
the tunneling currents for the ND and the N-DDW junc- 
tions. We start from the expression 



I(t) 



E i (^)( C l(^)C J , a (ky)) + h.C. (A 1) 



In Keldysh's formulation of non-equilibrium Green's 
functions, the expectation values (cj^c,,^) in the above 
equation can be expressed as a perturbation series. Un- 
der certain approximations one can resum this series to 
all orders in the tunneling amplitude t. 



1. ND junctions 

We consider first ND tunneling. In dealing with su- 
perconducting phenomena it is convenient to use the 
Nambu representation which explicitly distinguishes par- 
ticles and holes by assigning them to different compo- 
nents. One thus defines the spinor field-operators 



c a1; (xi,ky,t) (A 2) 



where a = {L, R} labels the electrodes and the upper 
and lower elements are associated with, respectively, elec- 
trons and holes. The Kcldysh_non-equilibrium Green's 
functions are then defined asHEa 

g^M^xj^') = + l {^lA x ^ t ')^*A^ t )) > ( A 3) 
5^(^^;^,0 = -*(*a, M (^,t)*^(^,t')) . (A 4) 

For brevity we have suppressed the k y dependence. The 
Green's functions here carry the left right indices a, (3 — 
{L, R}, the Nambu (spinor) indices [i,v = {1, 2}, and the 
Keldysh indices {— , +}. For notational clarity we shall in 
the following frequently omit irrelevant indices and keep 
track of only those related to our discussion. 

In this representation we define the tunneling matrix 



t = tr 3 a 3 



(A 5) 



where t 3 and a 3 are the third Pauli matrices pertaining 
to the Nambu space and the Keldysh space, respectively. 
In particular, 173 is chosen so that in the Keldysh space 



1 



, and (To 







since we have assigned the forward time-path the "— " 
time axis, and the return time-path the "+" time axis. 
In the following we will consider only real valued t and 
hence t* = t. 



The current expression (A 1) can now be written as 
J-oc 27r 

x { TT [g R l(xo,ky,Lu)} - TT[gl+(x ,k y ,u)]j . (A 6) 

where the trace is taken over the Nambu space and we 
have Fourier transformed over the time variables. In the 
presence of particle-hole symmetry the Nambu compo- 
nents 11 and 22 in the trace above contribute equally. 
Therefore, the trace yields twice the contribution— from 
the 11 component. Applying the Dyson cquationsEil 



go + go t g = go + g t g 



(A 7) 



and noting that g RL = = g LR , as there is no propa- 
gation between the two electrodes at the bare level, one 
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finds 



9rl = [9o + gigo] R t 
[go + go i 



*(SflflSo,Li " 9rr9o,1 l ) (A 8) 

•Jiii : [9a+gatg\ LR 

K-9o,tL9RR + 9o,Il9rr) • (A 9) 



Since the components of Keldysh Green's functions has 

- g H . Taking the diffcr- 



9 



the property g + g~ 

ence between g R £ and g^ R iising Eqs. (A 8) and (A £) 
and inserting the result into (A 6) we obtain 



9oXlm^ ~ eV )9RR,n{u) 

-gt,lL,u(P ~ eV)ff^ u (w)} -(A 10) 

Note that the frequency arguments of the bare Green's 
functions for the left electrode g has been shifted 

due to the applied bias eV between the two sides (/i^ — 
Hr = eV). We emphasize that at this stage Eq. ( A 1C ) 
is exact and our remaining task is to express the ex- 
act Green's functions g~RR^n + in terms of the uncoupled 
Green's functions go via some approximations. 

The Green's functions g^ / ^ can be expressed in 
terms of the bare Green's functions g$ ^ + and the 
exact Green's functions g r ' a by means of the following 
equationai-1 



9 



7- 



+-/- 



[tVM + 1] • (A ii) 



To fully reduce to bare quantities, one can express the 
exact Green's functions g r,a in terms of the bare Green's 
functions <?g' a by virtue of the Dyson equations 



(A 12) 



Ignoring many-body effects, we are able to resum the 
series and obtain 



9rr(u) = T RL (w) go,R R {u) , 
ffz&M = T lr^) \t gl[ a LL (uJ - eV)T 3 g r Q - a RR (uj) 
= Trl(u) [t g^ RR {uj)T 3 g^ LL {uj - eV) 
9ll{u) = TLr{u) 5o'1l(w - eV) , 



(A 13) 
(A 14) 

(A 15) 
(A 16) 



where sum over tunneling processes of all orders is signi- 
fied by the factors 



TrlM = 1 ~ t 9o,rr( u ) t 39o%(^ ~ eV)r 3 



q-r,a 
LR 



1 - 1 VoIl^ ~ eV ) T 39o% R (^)T3 



Since /zl — (1r — eV all frequency arguments of go.LL are 
to be taken at lu — eV. In arriving at the above equations 
we have used for the normal electrode 



go,LL{u-eV) = ( o 





50,LL,22(^ + eV) 



Note that the 11 and the 22 elements of <?o,ll represent 
particles and holes, respectively, and hence their response 
to applied bias is opposite. This is essential in giving rise 
to the Andreev contrib utions i n the t unnelin g curr ent. 
Incorporating Eqs. flA 13[ )-( [A 16| ) with ( A ll[ ), one 

can thus obtain g R ^ + and substitute back into (A 10). 
Finally we make use of the following relations 



«? -+H = 2«/( w )i( w ), (A 17) 

g+~(u) = -2Tn[l-f(u)]A(u), (A 18) 

where f(ui) is the Fermi function and A(uj) is the spectral 
weight matrix given by ( j3l"l) . This leads us to the current 
expressions (p7|)~(p0|). 



2. N-DDW junctions 

We turn now to deriving the current expressions for N- 
DDW junctions which is also applicable to NG junctions. 
We shall also show that in this case Andreev-like pro- 
cesses do not contribute to the tunneling current. In the 
absence of external fields, spin degree of freedom merely 
introduces a factor of two. Thus the spin indices a will 
be omitted in the following. 

We fi rst define the Keldysh Green's functions similarly 
to flO) 



9 a p{xi,t;xj,t') = +i(c^{x j ,t')c a {x u t)\ , (A 19) 
9ip{Xi,t;xj,t!) = -i(^ a (Xi,t)c^(Xj,1?)\ . (A 20) 

Here the subscripts a, j3 = {R, L} are labels for the elec- 
trodes (not to be confused with the labels for sublattices 
in the text). In terms of t he Keldysh Green's functions 
the tunneling current (|A l| ) can be written 



1(f) = +e <£+(r, I) - t* 9lt{l, r)] . (A 21) 



l.r,a 



Similar to the previous section, the renormalized Green's 
functions g RL and g LR can be expressed as combinations 
of the bare Green's functions g LL and the renormalized 
Green's function g RR . This results in the exact formula 

1 = e E ^ i £ Uo,ll (<•> - eV)g+ R (u;) 

Ky ,fT 

-g+l L (io - eV)g R +(io)} . (A 22) 

Note that here, unlike the previous section, the tunneling 
matrix is 



t = t<T3 . 



(A 23) 



where 173 is the third Pauli matrix in the Keldysh space. 

To proceed, as in the ND case, we apply ( [A 11| ) and 
expand the exact Green's functions g^R + in terms of 
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the bare Green's functions g$ ^ + and the renormalized 
retarded and advanced Green's functions g r ' a . Again, 
ignoring many-body effects we can express the renormal- 
ized Green's functions g r - a in terms of th e bare Green's 
functions g r Q ,a . Utilizing ( A 17 ) and ( A IE ), one has 



g R + = 2m[f(u)M R (u) + f(u> - eV)M L (u)] 

g+ n = -2m[(l ~ f(u))M R (u) + (1 - f(u - eV))M L (u)\ 

with 

M r (uj) = A R (u)\l+tg RL (uj)\ 2 , (A 24) 

M L {u) = t 2 A L (u;-eV)\g RR (u;)\ 2 . (A 25) 

In the last expressions we have used g r RL = {g% R )* ■ Note 
that Ml contains the spectral weight Al of the normal 
electrode. It is associated with tunneling processes where 
a particle is reflected back into the left side and at the 



same time a particle-hole pair is transmitted into the 
right; this is reminiscent of the Andreev channel in ND 
tunneling (cf. the integrand in (Boh). Su bstitu ting the 
above results into the current expression ( A 22 ), we ob- 
tain for the terms in the braces in the integrand 

^ 2 [f{uo - eV) - f{u)]A L {u ~ eV)M R (u) , (A 26) 



which leads to the current expression (|4£ 

It is remarkable that Ml is canceled completely in the 
final current expression, and hence no contribution from 
the Andreev process remains in the tunneling current. 
The reason for this difference between the ND and the N- 
DDW (and likewise NG) tunneling is that for the former 
the bias voltage shifts the energy of particles and holes 
in opposite directions, while for the latter particles from 
different bands acquire the same shift under applied bias. 



1 I. Giaever, Phys. Rev. Lett. 5, 147, 464(1960). 

2 E. L. Wolf, Principles of Electron Tunneling Spectroscopy 
(Oxford University Press, New York, 1985). 

3 C-R Hu, Phys. Rev. Lett. 72, 1526 (1994). 

4 This approach was pioneered by G. E. Blonder, M. Tin- 
kham, and T. M. Klapwijk in Phys. Rev. B 25, 4515 (1982), 
where they considered tunneling into s-wave superconduc- 
tors. The d-wave version was accomplished by S. Kashi- 
waya, Y. Tanaka, M. Koyanagi, H. Takashima, and K. Ka- 
jimura, ibid 51, 1350 (1995). 

5 See for example Y. Tamura et al., Phys. Rev. B 60, 9817 
(1999). 

6 M. P. Samanta and S. Datta, Phys. Rev. B 57, 10972 
(1998). 

7 C. L. Wu, PhD thesis, National Tsing Hua Univ (Hsinshu, 
Taiwan, 2000), and C. L. Wu, C.-Y. Mou, and D. Chang, 
Phys. Rev. B 63, 172503 (2001). 

8 J. C. Cuevas, A. Martfn-Rodero, and A. Levy Yeyati, Phys. 
Rev. B 54, 7366 (1996). 

9 The Green's function approach has also been applied at 
the quasi-classical level; see, for example, T. Luke et al., 
Phys. Rev. B 63 064510 (2001). 

10 C. Caroli, R. Combescot, P. Nozieres, and D. Saint- James, 
J. Phys. C: Solid St. Phys. 4, 916 (1971). 

11 X. Z. Yan, H. W. Zhao, and C. R. Hu, Phys. Rev. B 61, 
14759 (2000). 

12 The dimensionality of 2 is not essential; it is chosen here 
mainly for current experimental and theoretical interests. 
Our formulation is equally applicable to 3D cases. 

13 G. D. Mahan, Many Particle Physics 2nd ed. (Plenum, New 
York, 1990) p788. 

14 L. V. Keldysh, Sov. Phys. JETP 20, 1018 (1965). 

15 See, for example, J. D. Jackson, Classical Electrodynamics 



2nd ed. (John Wiley & Sons, N ew York, 1975) 



S-T Wu and C.-Y. Mou, arXiv: cond-mat/0112262 



Here, for simplicity, w e ha ve omitted the labels for elec- 
trodes {L, R} (cf. Eq. (A 2) in the Appendix). 
Y. S. Barush and A. A. Svidzinsky, Phys. Rev. B 55, 15282 
(1997). 

M. Fogelstrom, D. Rainer, and J. A. Sauls, Phys. Rev. 
Lett. 79, 281 (1997). 

P. G. de Gennes, Superconductivity of Metals and Alloys 
(Addison- Wesley, New York, 1989). 

M. Aprili, E. Badica, and L. H. Greene, Phy. Rev. Lett. 
83, 4630 (1999). 

Y. Dagan and G. Deutscher, Phy. Rev. Lett. 87, 177004 
(2001). 

W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. 
Lett, 42 1698 (1979). 

In reality the bond length in polyacetylene alternates with 
a variation of about 6% in the single and the double bonds. 
See, for example, T. Timusk and B. Statt, Rep. Prog. Phys. 
62 61 (1999) and references therein. 

S. Chakravarty, R. B. Laughlin, D. K. Morr, and C. Nayak, 
Phys. Rev. B 63 094503 (2001). 

I. Affleck and J. B. Marston, Phys. Re v. B 37 3774 (1988) . 
C. Honerkamp and M. Sigrist, arXiv: |sond-mat/010755C . 
M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, 
J. Phys. Soc. Jpn. 65, 1920 (1996). 
S-T Wu and C.-Y. Mou, in preparation (2002). 
X. Z. Yan and H. Iyetomi, Phys. Rev. B 57, 7944 (1998). 
J. B. Ketterson and S. N. Song, Superconductivity (Cam- 
bridge University Press, Cambridge, 1999) pp 380-383. 
A. Levy Yeyati, A. Martm-Rodero, and F. J. Garcia- Vidal, 
Phys. Rev. B51, 3743 (1995). 



